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Abstract. Within the framework of Von Neumann's expanding model, we study the 
maximum growth rate p* achievable by an autocatalytic reaction network in which 
reactions involve a finite (fixed or fluctuating) number D of reagents, p* is calculated 
numerically using a variant of the Minover algorithm, and analytically via the cavity 
method for disordered systems. As the ratio between the number of reactions and that 
of reagents increases the system passes from a contracting (p* < 1) to an expanding 
regime (p* > 1). These results extend the scenario derived in the fully connected 
model (D ~ > oo), with the important difference that, generically, larger growth rates 
are achievable in the expanding phase for finite D and in more diluted networks. 
Moreover, the range of attainable values of p* shrinks as the connectivity increases. 



1. Introduction 

Von Neumann's expansion problem was initially formulated to describe growth in 
production economies as an autocatalytic process and has played a key role in the 
development of the theory of economic growth [1-3]. In a nutshell, it concerns the 
calculation of the maximum growth rate achievable by an autocatalytic system of M 
reactants (labeled by Greek indices like /i, z/, . . .) interconnected by iV chemical reactions 
(labeled by Roman indices like i, j, . . .) specified by a given stoichiometric matrix. The 
basic ingredients are however sufficiently simple to be applicable in different contexts 
ranging from economics to systems biology. 

In order to state the optimization problem in mathematical terms (see however [4] 
for a more detailed description), it is convenient to separate the matrix of input 
stoichiometric coefficients A = {af > 0} (input matrix for brevity) from that of outputs 
B = {6f > 0}. The relevant microscopic variables are the reaction fluxes, denoted by 
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Si, which are assumed to be non- negative. The problem amounts to finding a vector 
s = { s i > 0} of positive fluxes and a number p > such that p is maximum subject to 



where the inequality is to be understood component-wise i.e. valid for each reactant p. 
The trivial null solution s = 0, viz. {sj = Vi}, is not accepted: at least one component 
of s must be strictly positive. 

Condition (1) requires that the total output of every reactant is at least p times 
the total input. So the maximum feasible p, p*, measures the largest uniform growth 
rate achievable by the given set of reactions. If p* < 1, the optimal state of the system 
is a contracting one. Otherwise for p* > 1 it is expanding. The case p* — 1 describes 
instead a system in which at optimality reaction fluxes are arranged so as to guarantee 
mass balance. 

Von Neumann's problem has been studied recently in a statistical mechanics 
perspective under the assumptions that M scales linearly with iV (with iV — > oo) and 
that the stoichiometric matrices A and B have quenched random independent and 
identically distributed elements [4]. With a fully connected network of reactions (that 
is, one in which every reaction uses each reactant both as an input and as an output), 
one finds a transition from a contracting to an expanding regime when the parameter 
n = N/M exceeds the critical value n c — 1. 

In this work we extend the analysis of [4] to the finitely connected case where 
reactions use a finite number of inputs to produce a finite number of outputs. By 
analogy with integer programming problems (see e.g. [5]), we represent our autocatalytic 
network as a bipartite (factor) graph (see Fig. 1). We denote reagents as squares and 
reactions as circles. Each circle has a number of incoming (resp. outgoing) connections 
to squares, representing the inputs (resp. outputs) of the reaction, and each of these 
links carries the input (resp. output) stoichiometric coefficient. Similarly, each square 
has a number of incoming (resp. outgoing) connections to circles, denoting the reactions 
where that reagent enters as an output (resp. input). The notation % e p identifies a 
reaction % in which a reactant p is involved either as an input or as an output (and 
vice- versa for p e i). To each reaction node a variable > is attached, representing 
the flux of reaction i. Each reactant node p carries instead a function c M given by 



We are interested in finding non-null flux configurations satisfying (1), that is such 
that & 1 > for each p and specifically in finding the largest value of p for which a 
configuration of this type exists. Our approach is developed along two main lines. 
On the one hand (Section 2), we employ a suitably modified Minover algorithm to 
compute optimal growth rates numerically on any graph. Results thus derived will be 
theoretically validated in Section 3, where we use the cavity method to describe solutions 
analytically for locally tree-like instances. This approximation turns out to be in good 



Bs > pAs 
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Figure 1. Factor graph representation of a network of reactions and reagents. 
Reactions are represented by circles, reagents by squares. Stoichiometric coefficients 
label the directed links. 

qualitative agreement with numerical data. The connection of the latter approach with 
the fully connected theory developed in [4] is presented in two detailed Appendices. 

2. Numerical analysis 

2.1. Algorithm 

Formula (1) is reminiscent of pattern storage conditions in the perceptron, a formal 
neural network model [6] . This analogy allowed us to devise an exact algorithm to solve 
Von Neumann's problem on any graph in a time growing only polynomially with the 
number N of nodes. The algorithm is an extension of the so-called Minover algorithm [7] 
enforcing the constraint of positivity for the fluxes, and hereafter referred to as Minover + . 

The core part of the algorithm is a subroutine, Minover + (p), we now describe. We 
assume that the matrices A, B are given, and p denotes a real, non-negative number. 

(i) Step £ = 0: all fluxes are null, Sj = 0. 

(ii) Step £ + 1: Calculate the M functions c M (£) defined in (2) from the values Si(£) of 
the fluxes at step £, and look for the reactant with least value, 

Ho(t) = arg minc^) ; (3) 

* if c"°M > then OUTPUT s{£) and HALT (the null solution is not accepted); 

* otherwise, if c w( ^ < 0, update the fluxes according to the rule 

s t {£ + 1) = max{0, s t (£) + 6f oW - paf oW } , (4) 
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increase £ by one, and GO TO (ii). 

In case of ties, one may select p at random uniformly among those with the same 
(and lowest) values of & '. 

We claim that, if p < p*, the largest growth rate associated to matrices A and B, 
Minover + (p) will halt after a finite number of steps, and output a set of non-negative 
fluxes guaranteeing a growth rate larger than p. The proof is simple and goes as follows. 
Assume p < p* . Then there exists a strictly positive number A, hereafter called stability, 
and a vector s* ^ of fluxes such thatf 

s* • (fe" - pa") > A |s*| for each p . (5) 

Let 

A = max5>f-paf) 2 >0. (6) 

% 

and consider the functions M{£) = s* ■ s(£) and T>(£) = \s(£)\ 2 . We have, from (4), 
M(t + 1)=5>* max{0, Si (£) + 6f o(£) - paf oW } 

i 

> E^^w+ fe " oW -p«r w ) (7) 

i 

> Af(£) + A \s*\ 

from the non- negativity of s* and (5). We thus obtain M(£) > £ A \s*\ for any step £. 
Turning to function V(£) we get 

V(£+l) = ^max{0, Sl (£) -paf oW } 2 

% 

< E(^w+ & r w -P«r w ) 2 (8) 

i 

= V{£) + 2c«'W + ^(6f oW " P<° {i) ? 

i 

< V(£) + A 

from definition (6) and the assumption that the algorithm has not halted at step £ i.e. 
c»W < 0. Hence V(£) < £ A. Now let 



By Cauchy-Schwarz inequality f(£) < 1. But according to the above calculations 

m = -tm => ^ =yrt . do) 

Therefore the algorithm halts after 

^o = ^ (11) 

| The set S(p) of fluxes fulfilling constraints (1) is, by definition of p* , non empty. Any vector s* in 
the interior of S(p) satisfies (5) for some positive A. A natural choice is the optimal vector of fluxes, 
i.e. that associated to p*. 
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Figure 2. Convergence time (in individual steps) versus p for a single network with 
M = 100 reagents and n — 0.66. The number of reagents per reaction is fixed at 
D = 5, the number of reactions per reagent is a Poisson random variable with mean 
Dn = 3.3. The dashed line represents the best fit with r conv ~ (p* — p)~ 2 . For this 
network, p* = 0.5994. 



steps at most. Now call £ s the step at which the algorithm stops. The halting condition 
ensures that c^ ^ > 0, i.e. condition (1) is satisfied, with a non-null vector of fluxes 
(by construction the fluxes are non-negative at any step). Hence s(£ s ) is a solution to 
Von Neumann's problem with growth rate p. 

Our complete algorithm, Minover + , is an iteration of Minover + (p) for increasing 
values of p. Starting from p = e, a small positive value for which (1) surely admits a 
solution §, one can measure the convergence time (number of steps prior to halting) £ s for 
increasing values of p. The convergence time increases too. Now we expect A to vanish 
as p* — p when p gets closer and closer to its optimal value, thus £ s ~ (p* — p) -2 from 
(11) (see Fig. 2). One may then estimate the maximal growth rate by extrapolating 
from the log-log plot of the convergence time versus p. 



2. 2. Survey of results 

In what follows we restrict ourselves to purely autocatalytic systems, that is we assume 
that every reagent is produced and consumed by at least one reaction. Indeed a reagent 
that only serves as input gives rise to a constraint of the form 

c M = -p^^af>0 (12) 

which is satisfied by taking p = 0. Such a situation would immediately force the result 
p* — 0||. On the other hand, "sink" reagents (namely those that are not inputs of any 

§ We can always choose e = min i!M 6f/af > 0. 

|| The a priori probability to generate a factor graph without such pathologic function nodes in the 
Poissonian model discussed later is given by (1— e~ Dn (\~ e~ Dn )) M , where M is the number of reactants. 
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Figure 3. Left panel: p* vs n for Regular(iV)/Poisson(Af) networks of 100 reagents 
for various D. The continuous line is the fully-connected limit for a system of the same 
size [4] while F.C. labels numerical results for a fully connected system. Right panel: 
p* vs n for Regular (iV)/Poisson(M) (open markers) and Regular (N) /SF(M) (closed 
markers) networks for different values of D (see text for details). Averages over 50 
samples. 



reaction) provide constraints that are p-independent and trivially satisfied. Furthermore, 
for the sake of simplicity, we fix the number M of reagents and assume that the quenched 
random stoichiometric matrix has Gaussian elements with mean 1 and variance 1/2. 

For a start, we consider Von Neumann's problem on factor graphs in which every 
reaction has D inputs and D outputs, while for reagents we assume that the in(out)- 
degree distribution is Poissonian, P(k) = X k e~ x /k\ with mean degree A = Dn. For 
short, we write Regular(iV)/Poisson(M) to denote this type of situation. Our analysis 
focuses on values of n = N/M > 1/D, which ensure that the resulting factor graph is 
connected. Results for p* in this case are displayed in the left panel of Fig. 3. One 
sees that the overall qualitative behaviour reproduces the results obtained for the fully 
connected model. The system passes from an expanding (p* > 1) to a contracting 
(p* < 1) phase when n is lowered below a critical value. Interestingly, the maximum 
growth rate achievable in the diluted system is larger than in the fully connected model 
in the expanding phase and higher growth rates can be achieved in more diluted systems. 
This conclusion turns out to be valid generically for all types of networks we studied. 

In the right panel of Fig. 3 we compare these results with those obtained for the case 
in which the number of reactions per reagent is power-law distributed (SF or scale-free 
for short) rather than Poisson, which enables the coexistence of widely used and rarely 
used reagents in the reaction network (reactions still have a ^-distributed connectivity). 
More precisely, the number of reactions per reagent is distributed as P(k) ~ k" 1 with 
2 < 7 < 3 (as in e.g. metabolic networks [9]). The particular value of the exponent does 

This implies that such nodes are always present in the thermodynamic limit. In this work we discard 
them completely. However in practical applications (e.g. metabolic networks [8]) it is important to 
take these nodes into account as prescribed "sources" (uptakes). 
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Figure 4. Left panel: p* vs n for Poisson(iV)/Poisson(M) networks. Right panel: p* 
vs n for different types of factor graphs with D = 5. Averages over 50 samples; system 
of 100 reagents. 



not affect results in a marked way. Note that a significant difference (which however 
lies inside the error bar) is observed only for the smallest values of D. 

Next we compare networks with degree regular reactions, where, as said above, the 
in/out-degree of reactions is fixed equal to D, with Poisson(iV)/Poisson(M) networks 
where also reactions have fluctuating degree with mean D. From Fig. 4 one sees 
that optimal growth rates for all values of n are always larger in Poissonian networks 
than in regular ones, independently of the degree distribution of reagents (right panel). 
Moreover a fluctuating connectivity for reactions makes the expanding regime achievable 
with fewer reactions (left panel), and in particular the critical point where p* becomes 
larger than 1 is n c ~ 1, similar to what is found in the fully connected case. In few words, 
we can say that topological regularity has a strong influence on the maximum growth 
rate achievable and that an asymmetric choice of degree distributions for reactions and 
reagents moves n c from its fully connected value. In the case above, it takes a larger 
repertoire of reactions to sustain an expanding regime. 

It is possible to find a simple expression that allows to re-scale the data obtained 
for different -D's onto a single curve (which appears to be topology-dependent), at least 
for large n. It turns out that 



where pp C (n) is the optimal growth rate of the fully connected system (see Fig. 5). 
This scaling form, which gives p*'s for diluted systems as 1/D corrections to the fully 
connected limit, implies that the range of achievable growth rates shrinks as D increases. 

In Fig. 6, we show the values of p* for systems of different sizes and topologies 
(similar results hold for other types of networks). One sees that already for moderate 
system sizes finite size effect are negligible. 

To conclude, we present (see Figures 7 and 8) the flux distribution at optimality for 
Regular(iV)/Poisson(M) and Poisson(iV)/Poisson(M) networks for values of n below 
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Figure 5. Scaling functions /(n) for Regular(7V)/Poisson(Af) networks (left) and 
Poisson(iV)/Poisson(M) networks (right). Error bars not reported. 
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Figure 6. p* vs n for Regular(A r )/Poisson(M) and Regular(iV)/SF(M) networks of 
100 and 200 reagents for different D. Averages over 50 samples. 



and above the critical point. In the former case, distributions are well fitted by 
an exponential form both for small and large n. In the latter case, an algebraic law 
provides the best fit in the expanding regime. Similar results are obtained for different 
connectivity distributions for reaction nodes. The particular value of the exponent 
appears to be topology dependent. However, it is interesting that such distributions arise 
in this context, specifically for their relevance in connection to metabolic networks [8,9]. 
Note finally that as in the fully-connected model the weight of the mass at s = 
increases with n, though less drastically 



P(s) 



0.1 



0.01 



0.001 



0.0001 



I I I 1 1 1 II 1 — I I I I 1 1 1 1 1 — I I I 1 1 1 II 



n=0.66 



\ 



(|) = 0.213 



V 

m 
•a 



I i I i i i i lME : 

0.1 1 10 

s 



j n if 1 — i i i 1 1 in 1 — i i i i 1 1 1 1 1 — i r 



r 

n=3.33 



0.1 - 



0.01 E~ 



0.001 



0.0001 



4 = 0.257 



\ 



1 1 1 1 1 i i i 



J m ■ i ■ i 



0.1 



10 



Figure 7. Flux distribution for a Regular(A r )/Poisson(M) network of reagents with 
D = 5 inputs and outputs per reaction and different values of n (average over 50 
samples). System of 200 reagents. The 5-peak at s = (carrying an intensive weight 
<p reported inside the panels) is not shown. 




Figure 8. Flux distribution for a Poisson(7V)/Poisson(Af ) network of reagents with 
D = 5 inputs and outputs per reaction and different values of n (average over 50 
samples). System of 200 reagents. The <5-peak at s — (carrying an intensive weight 
</> reported inside the panels) is not shown. The red lines represent best fits with an 
exponential distribution of the form p(x) — Ae~ Bx (left panel) and with an algebraic 
distribution of the form p(x) = A(B + x)~ c (right panel). 



3. Cavity theory 

Von Neumann's problem possesses the natural cost function 

M 

that, given p, counts the number of reagents for which the condition (1) is violated. It 
is clear that p* corresponds to the largest p for which one can arrange fluxes in such a 
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Figure 9. Left: reagents participating in reaction i are correlated. Right: If reaction 
i is removed those reagents become uncorrelated on a tree-like graph. 



way that E = 0. More generally, we can look for configurations that minimize E, so it 
is convenient to introduce the "partition function" 

Z = J 1 (s)e- f3E ds = Tr s e~ pE (15) 

with 7(s) some constraint over the configuration space (e.g. a linear constraint). Minima 
of the cost function are selected in the limit (3 — > oo. 

Let us assume that the factor graph has a locally tree-like structure and consider, 
as in Figure 9, a reaction i and the reagents involved in it. It is clear that the joint 
distribution of consumptions c M£ * = {c^}^ does not factorize into single-reagent terms, 
i.e. 

q(c^)^l[ q »(cn (16) 

because the /x's share the common vertex i. However, if we remove the reaction % (here 
and in what follows indexes enclosed in parentheses like [x] denote quantities calculated 
removing the node x, while sums and products where x is not considered carry the index 
\x) the resulting joint distribution of consumptions q , (i)(c /ie? ) does factorize: 

Similarly, with Si £fl = {sj} igM , we have 

P(si^)^l[pi(si) (18) 

but, removing reagent 

P { "\s ie ,) = l[p^(s i ) (19) 
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Figure 10. From the graph where reaction i has been removed (left), one calculates 
easily the distribution p^- (si) in the absence of the reagent /i by merging back all 
reagents participating in reaction i but reagent fi. 



The quantities (?^(c M ) and (sj) are known as cavity distributions. Now for all 
% = 1, . . . , N and fi 6 % it is possible to find equations for the cavity distributions 
by simply merging back all disconnected nodes but one (see Fig. 10). It is easily seen 
in particular that 



1_ L-^ A ^\-^-sm- P <)\ -q ^)dc v 



1 



Tr. 



(20) 



ru 

v£i\fi 

(where is a normalization factor ensuring that J p^f\si)dsi = 1) while for each 
[A — 1, . . . , M and i e we have 



Tr, 



J 6 A* 



Note that these in turn imply that 



ft(«) 4^. 9(l )(^)e-^- e K)-^-K)] 
A: 



g^) = Tr s . >)( Sje/i )<5 



(21) 

(22) 
(23) 



so that the whole probability distribution can be reconstructed from the set of cavity 
distributions. The statistical average cost function is finally given by 

g"(c) [e-P + (1 - e-Q6 (c)] 
-c)— — r~/ \ r — s . a-, ~. / \ i (24) 



Tr c g^(c) [e-/ 3 + (1 - e~P)Q (c)] 
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The cavity theory developed here recovers the replica theory of [4] in the fully-connected 
limit. This is shown in detail in Appendices A, where a slightly revised replica theory 
is discussed, and B, where the fully connected limit is constructed explicitly. 

Because the configurational variables Sj and d 1 are continuous, solving equations 
(20) and (21), which in the zero temperature {(3 — > oo) limit take the form 



II e ft + *(&?- K)] (25) 
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(26) 



by the standard method of population dynamics requires studying a population of 
populations, as e.g. in [10] (in other words, it is equivalent to a one-step RSB calculation 
already at the replica-symmetric level). We therefore decided to focus our analysis on 
the calculation of the optimal growth rate, for which it is possible to obtain good results 
at a modest computational cost. 

We are thus interested in finding the critical line p*(n). One may proceed as follows. 
Assume there is no linear constraint (the value of p* is unaffected by the presence of a 
linear constraint). Then the equations always admit the trivial solution (s = 0) for all 
values of p. In fact, for p > p* this is the only zero-energy solution while below p* other 
non-null solutions exist. Thus, starting from random initial conditions for the fluxes, 
we would expect the average flux to vanish (resp. remain different from zero) under the 
iteration of the cavity equations for values of p above (resp. below) p*. This obviously 
assumes that the iteration of the cavity equations for p < p* keeps the fluxes away from 
their trivial value, which turns out to be the case. 

To check the validity of this assumption, we have first considered the cavity 
equations (25) and (26) in their fully connected limit (see Appendices for details), viz. 
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(27) 
(28) 
(29) 
(30) 

(31) 
(32) 
(33) 



C(p) = K - p< 

and where M s (a,b) and fk(a,b) are defined as follows (H(x) = erfc(x/v / 2)/2) 

(1-^) (D^M) 



M s (a, 6) = 



s=0 



fk{a,b) 



e-P + (1 - e-/ 3 

as—bs 2 c .k 



(/?->oo) 



Tr s e 
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(34) 
(35) 

(36) 
(37) 



Tr s e as ~ 6s2 

with £/ c (a, 6) a Gaussian distribution for the variable c with mean a and variance b. 
For computational reasons, it is more convenient to iterate the corresponding TAP 
equations [11], namely 



= M n 



m = fi 
h 



U 

B, 



(38) 
(39) 
(40) 

(41) 
(42) 
(43) 



Given a matrix £f(p), we solve the preceding set of coupled equations by iteration: 
starting from a small value of p we monitor the average flux and determine p* when, 
under iteration of the equations, the average flux goes to zero. In Fig. 11 we compare the 
results for a system of M — 100 reagents (average over 10 samples) with the theoretical 
prediction. The agreement is fairly good. 

We have then considered the same approach to calculate p* for a 
Regular(iV)/Poisson(M) network with connectivity D = 7 (see Fig. 11). To do so, 
we have considered the ensemble version of equations (25) and (26), which take the 
form (we suppress cavity indexes (x) for simplicity) 

„ 2D-1 



<f(c) 



•fcu-1 



Tr c , gV)e[c" + a (&r-K)] 



n ^'iPAsj 



kn-l 



(44) 
(45) 
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Figure 11. Results from cavity theory. The continuous (blue) line is the theoretical 
(replica) prediction in the fully connected limit. The dashed line represents Minover 
results for Regular(A r )/Poisson (M) graphs, with error bars. Markers correspond to 
iteration of the fully connected (circles, average over 10 samples, number of reagents 
M = 100) and Regular(iV)/Poisson (M) cavity equations (squares). 



with k p a random Poisson number with average nD. In the ensemble, which is obtained 
after the limit iV — > oo has been performed and where n plays the role of an external 
parameter, the indexes % and \l now take values i, p, — 1, 2, 3, . . . and describe an infinite 
population of probabilities {q^{c)} ^ = i^... and {Pi(s)}i=i,2,...- We have studied a finite 
population of probabilities, i.e. i,p = 1, . . . ,J\f pop , which are initialised randomly. In 
order to find p*, we start from a small initial value for p. The population is iterated 
via equations (44) and (45). If after a certain (sufficiently large) number of iteration 
steps the average flux is different from zero, we increase the value of p by a small 
value and restart the procedure. Eventually we can locate a value of p at which the 
average flux vanishes. This defines p*. In Fig. 11 we display the resulting curve p*(n) 
for a population of Af pop = 50 distributions and D = 7, which qualitatively agrees 
with the results found by the Minover + algorithm. (We remind that cavity equations 
are obtained for a tree-like topology and the network we considered in the numerical 
solution is relatively small.) 

4. Summary 

Despite evident quantitative differences, the overall picture obtained for the fully- 
connected Von Neumann's model is roughly preserved in the case of finite connectivity, 
discussed here. The major difference appears to be that much larger growth rates can 
be achieved in diluted networks (reactions with few inputs/outputs are evidently more 
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efficient) and, generically, when reactions have stochastic connectivities. The framework 
of Von Neumann's problem, including its global optimization aspect, is sufficiently 
simple to ensure that the methods developed in this work can be applied in several 
different contexts, particularly in economics and biology. Work along both lines is in 
progress. 
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Appendix A. The fully-connected model revisited 

We shall briefly re-consider here the replica approach of [4] for the fully connected model 
in the Derrida- Gardner framework (as opposed to the Gardner one employed in [4]). This 
allows us to derive equations for the order parameters that are easily compared with 
those which shall be derived in the following section as the fully connected limit of the 
cavity theory. 

Our starting point is the partition function, which we write in this case 
m m / n \ 

Z = Tr c J] [e-P + (1 - e-e)e{d>)] Tr s J] 5 I c» - -= £ s^{p) (A.l) 

fi=l /i=l \ V i=l / 

where 



= k - p< 



(A.2) 
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As in [4], we write fef = 6(1 + and af = a(l + af) with Gaussian af and /3f and 
further set a = 6 = 1. Moreover, we re-cast p as p = 1 + g/VN so that 

tf(p)=#-a?--^(l + a?) (A.3) 



The replicated and disorder-averaged partition function reads 

r M 

= iw.,c} n n + c 1 - e ^) e «)] 



= 1 /K=l 



r M / AT 



>< Tr {sl ,..., Sr} n n m - e ( a - 4 ) 



/ i /,- ! x v N i=l 



If the (^-distributions are written in their Fourier representation (with multipliers c^) the 
disorder average can be carried out easily, generating the standard order parameters 

qu ' = S Si£S ^' m ^ = a? S s ^ ( A,S ) 

i i 

which can be forced into the partition function with proper Lagrange multipliers qw 
and rh^. One finally finds 

*=f°«*™** ""^ < A6 » 

where (a = M/N, k = (/3f - a? ) 2 ) 



F(q,m) = a log ^ 



dcidfy 



n 



2tt 



x ;Q [e"' 3 + (1 - e-^e( Q )] 
S(q,m) = — logTr {sa e ~ E ^'«"' E -i s « s «'" E ^^ E »=i s « 



AT 

In the limit iV — > oo, the integral (A. 6) is evaluated by a saddle-point method. Let 
us consider the saddle-point equation for the order parameter q^ (I ^ £'): 

I He ^2? Wt' e i ^=^ ece+i9 ^ eme -^^'^' q "' nLi + (! - e^)0(c £ )] 

O//" — KCK ; r^rj -: — ^ ttzz (A. 

JYI £ dc ^ 5t e l T,e=i c e c e+ 1 9T,e^m e -^Y,e,e' c e c £'1W FJ^_ 1 [e~P + (1 — e _ ^)©(Q)] 

The argument of the exponential in the denominator can be linearized by a Hubbard- 
Stratonovich transformation. Once this is done, it is clear that the remaining integrals 
over q and q factorize and one easily sees that in the replica limit r — > the denominator 
tends to 1 and we can neglect it. It then remains to evaluate the denominator with the 
replica-symmetric (RS) Ansatz 

qu> = (Q - q)Sw +q m e = m (A. 8) 

qu> = (Q - q)5a> + q rhi> = m (A. 9) 
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in the replica limit r — > 0. With standard manipulations one finds 



q = —na I dxQ x (0,K,q) 



e 2k(Q- 9 ) » 



e-/ 3 + (1 - e^tf 



(A.10) 



y«co-9) / / 

where G x (a,b) denotes a Gaussian distribution for the variable x with mean a and 
variance 6. while H(x) = |erfc(a;/v / 2). 

Along similar lines one can derive the corresponding equations for the remaining 
order parameters: 

1 - e-P)d x Q x (-gm, k(Q - q)) 



Q 



dx G x (0, nq)- 



e~P + (1 - e~P)H f m+x 



m 



= S aJ 



dx Q x (0, nq)- 



_ (gm+x) 
g 2«(Q-q) 



"^ + (1 



-P)H 



gm+x 
n(Q—q) 



q= / dxQ x (0,l) F 

Q= dxg x (o,i)[± -—^ T 

f , ( fdse-(^> 2+{ix ^-^ s s 

m = / dxGJO, 1 i — — = 

J V fdse-( Q -* )' 2 +(**v^-* a )« 



(A.H) 

(A.12) 

(A.13) 
(A.14) 
(A.15) 



Appendix B. The fully-connected limit of the cavity theory for the diluted 
model 

In this Appendix we study the large connectivity limit of the cavity theory and derive 
equations for the relevant order parameters in this limit, in order to show that these 
equations are equivalent to those obtained in the revisited replica theory presented 
above. The first problem is to study the cavity distributions (20) and (21) in the fully 
connected limit for a fixed disorder realization. Next, we shall perform the disorder 
average. 

Let us begin from (21), which we re-cast as 



5g)(c) = Tr Sj6MV 



II r^te 



2tt 



(B.l) 
(B.2) 
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where we introduced a re-scaling factor 1/yK for the fields (K being the average 
connectivity). For K ^> 1 we can expand the exponential. Keeping terms up to the 
second order we find 



0) 



JJ Tr s .p5 M) ( Sj )e-^^ (p) = exp 
where 

m i M) = ( S i) 

^ = «> M - «, 

(/(*)>(,) = TW^M/M 
Inserting (B.3) into (B.2) and integrating over c one sees that for K 3> 1 



(B.3) 

(B.4) 
(B.5) 
(B.6) 



where as before G x (a, b) denotes a Gaussian distribution for the variable x with mean a 
and variance b. 

The limit of (20) is a bit more complicated. For a start, we re-write it by introducing 
the usual re-scaling: 



P?°(*) 



1 



n ^(o 



=^n^^( C )e-^h->^)] 



+ (i- e -/»)e(c + -=<(p) 



We focus on 

Tr.^^e-^h-^^^^Tr.^Cc) 
Now for any K the Heaviside function can be disposed of via 
Q(c+-±=s^(p))=Q(c)+e»(c,s) 



so that 



^ s ) = E^(^^r(p)) n ^- 1) ( C ) 

n>l ^ * ' 



(B.8) 
(B.9) 

(B.10) 

(B.H) 
(B.12) 



Tr cq ^(c)e-" e i- c -^ s ^} = Tr c ^(c) [e^ + (l - e^) 0(c)] 

+ (l-e-^Tr^f^cK^,,) 

Factoring out a term Tr c g^(c) [e _/3 + (l — e _/3 ) 6(c)] one sees that (B.9) becomes 

= i II [l + ll-e-^Tr.Q^^c,^] (B.13) 
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where is a new normalization factor and 



QUc) 



to 



(B.14) 



^ Tr e qfa(c)[e-fi + (l-e-fi)e(c)] 
With the expression given in (B.12), we have 

Tr c Q^(cK(c, S ) = fi^VirCjjj^lO) (B.15) 

For if 3> 1 we can truncate the above expansion after the first two terms and 
approximate the resulting expression with an exponential. This gives us 



A^s-B^s 2 



Tr, e 



(B.16) 



with 



= (1 - e^) ^ £ [£f (p)f (0) + (1 - e"') [Q^)(0)] 2 )(B.18) 

(B.7) and (B.16) represent the K ^> 1 limits of the cavity distributions for a fixed 
disorder sample. 

Recalling (A. 2) and (A. 3), we can now evaluate averages over the quenched disorder. 
To this aim, we set p = 1 + g/VK so that 



t 9 



(B.19) 



Computing the statistics of the quantities and over disorder one obtains 
(with a = K/N) 



Af } =A=-ga(l- e~ p ) F 



(4 M) - A) 2 = a A = aK{l- e-Py D 



= B = an (1 - e-P) [C + (l - e^) D] 



(B.20) 
(B.21) 
(B.22) 



(B^ is sample-independent in the fully connected limit K — > oo, so we neglect its 
fluctuations) where we used the following shorthands: 



c = 


<% (1) (o) 


D = 


[2w(o)] 2 


F = 




K = 


w - <y 



This implies that 



[</(*)>(„)]* 



Tr s /(s)e 



Hs-Bs 2 



(B.23) 

(B.24) 
(B.25) 
(B.26) 

(B.27) 
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We now must evaluate C, D and F, recalling Q^(c) as defined in Equations (B.14) and 
(B.7). The statistics of h = ^ E ie ^^(p)^ } is given by 



h = -gM, 

¥ = g 2 M 2 + nq, 



M = m 



Q = (si 



(h — h) 2 = Kq 
For the variance of Q v ,~ (c) we obtain instead 



K 

We can therefore write 



2 \ 

mm 



F 
D 

C 



J dhG h (-gM, nq) 
J dhg h (-gM, nq) 
J dhg h (-gM, nq) 



G c =o [h, k(Q - q)} 



J dc G c [h, k(Q - q)} {e-P + (1 - e-P) (c)] 

Qc=o [h, k{Q - q)} 

/ dcQ c [h, k(Q - q)} [e-P + (1 - e~P) (c) 

Gill, [h, k{Q - q)] 



J dc g c [h, k(Q — q)] [e~P + (1 - e-P) (c)] 
These form a closed system together with the equations for M, Q and q: 

Tr s s e Hs - Bs ' 2 



M = T^) = J dH g H (A,a A ) 
Q = W^) = J dH g H (A,a A ) 
q = ( s »)(ju) = J dH Gh(A, a a ) 



Tt s e 1 



Hs-Bs 2 

Tr s s 2 e Hs - Bs ' 
Ti s e Hs ~ Bs2 

Ti s s e Hs ~ Bs2 

Tl Hs-Bs* 



B.28) 
B.29) 
B.30) 

B.31) 

B.32) 
B.33) 
B.34) 

B.35) 
B.36) 
B.37) 



With minor manipulations, these equations can be precisely identified with the equations 
obtained by the replica approach in RS approximation in the previous section. 



